package model.util;

public class KirkJones {

	public enum Carrier {
		SINUSOIDAL, RECTANGLE;
	}

	public enum Order {
		DIRECT, INVERSE;
	}

	private double uSize;
	private double vSize;
	private double uFrequence;
	private double vFrequence;
	private Carrier carrier;
	private Order order;
	private double maximumValue;
	private double[][] matrixRe;
	private double[][] matrixIm;
	private double[][] hologramMatrix;

	public KirkJones(float uSize, float vSize, float uFrequence, float vFrequence, Carrier carrier, Order order,
			double[][] matrixRe, double[][] matrixIm) {
		this.uSize = uSize;
		this.vSize = vSize;
		this.uFrequence = uFrequence;
		this.vFrequence = vFrequence;
		this.carrier = carrier;
		this.order = order;
		this.matrixRe = matrixRe;
		this.matrixIm = matrixIm;
		hologramMatrix = new double[matrixRe.length][matrixRe[0].length];

	}

	public double[][] perform() {
		maximumValue = 0;

		for (int i = 0; i < matrixIm.length; i++) {
			for (int j = 0; j < matrixIm[0].length; j++) {
				double mod = Math.sqrt(matrixIm[i][j] * matrixIm[i][j] + matrixRe[i][j] * matrixRe[i][j]);
				maximumValue = maximumValue < mod ? mod : maximumValue;
			}
		}

		for (int i = 0; i < matrixIm.length; i++) {
			for (int j = 0; j < matrixIm[0].length; j++) {
				if (carrier.equals(Carrier.SINUSOIDAL) && order.equals(Order.DIRECT)) {
					double maximum = 2.38 + Math.PI * 2;
					double minimum = -2.38;
					double arg = 0;
					double mod = Math.sqrt(matrixIm[i][j] * matrixIm[i][j] + matrixRe[i][j] * matrixRe[i][j]);

					if (matrixIm[i][j] == 0) {
						arg = Math.PI / 2;
					} else {
						arg = Math.atan2(matrixRe[i][j], matrixIm[i][j]);
					}

					double u = (uSize / matrixRe.length) * (i - matrixRe.length / 2);

					double v = (vSize / matrixRe[0].length) * (j - matrixRe[0].length / 2);

					hologramMatrix[i][j] = arg + bessel0(mod / maximumValue)
							* Math.cos(2 * Math.PI * (uFrequence * u + vFrequence * v));

				}
				if (carrier.equals(Carrier.SINUSOIDAL) && order.equals(Order.INVERSE)) {
					double maximum = bessel1(0) + Math.PI * 2;
					double minimum = -bessel1(0);
					double arg = 0;
					double mod = Math.sqrt(matrixIm[i][j] * matrixIm[i][j] + matrixRe[i][j] * matrixRe[i][j]);

					if (matrixIm[i][j] == 0) {
						arg = Math.PI / 2;
					} else {
						arg = Math.atan2(matrixRe[i][j], matrixIm[i][j]);
					}

					double u = (uSize / matrixRe.length) * (i - matrixRe.length / 2);
					double v = (vSize / matrixRe[0].length) * (j - matrixRe[0].length / 2);

					hologramMatrix[i][j] = arg + bessel1(besj1(Math.PI / 2) * mod / maximumValue)
							* Math.cos(2 * Math.PI * (uFrequence * u + vFrequence * v));

				}

				if (carrier.equals(Carrier.RECTANGLE) && order.equals(Order.DIRECT)) {

					double maximum = Math.PI * 3;
					double minimum = -Math.PI;
					double arg = 0;
					double mod = Math.sqrt(matrixIm[i][j] * matrixIm[i][j] + matrixRe[i][j] * matrixRe[i][j]);

					if (matrixIm[i][j] == 0) {
						arg = Math.PI / 2;
					} else {
						arg = Math.atan2(matrixRe[i][j], matrixIm[i][j]);
					}

					double u = (uSize / matrixRe.length) * (i - matrixRe.length / 2);
					double v = (vSize / matrixRe[0].length) * (j - matrixRe[0].length / 2);

					hologramMatrix[i][j] = arg + Math.acos(mod / maximumValue)
							* Math.pow(-1, uFrequence * u + vFrequence * v);
				}
				if (carrier.equals(Carrier.RECTANGLE) && order.equals(Order.INVERSE)) {

					double maximum = Math.PI * 3;
					double minimum = -Math.PI;
					double arg = 0;
					double mod = Math.sqrt(matrixIm[i][j] * matrixIm[i][j] + matrixRe[i][j] * matrixRe[i][j]);

					if (matrixIm[i][j] == 0) {
						arg = Math.PI / 2;
					} else {
						arg = Math.atan2(matrixRe[i][j], matrixIm[i][j]);
					}

					double u = (uSize / matrixRe.length) * (i - matrixRe.length / 2);
					double v = (vSize / matrixRe[0].length) * (j - matrixRe[0].length / 2);

					hologramMatrix[i][j] = arg + Math.asin(mod / maximumValue)
							* Math.pow(-1, uFrequence * u + vFrequence * v);
				}
			}
		}
		return hologramMatrix;
	}

	double bessel0(double arg) {

		// ************ Transcendent bessel equation solution
		// **************************

		double arg1;

		arg1 = arg;

		if (arg > 1)
			arg1 = 1;
		if (arg < 0)
			arg1 = 0;

		double a[] =

		{ 2.38, 2.37, 2.355, 2.335, 2.32, 2.3, 2.285, 2.265, 2.25, 2.23, 2.21, 2.19, 2.175, 2.16,

		2.145, 2.13, 2.105, 2.085, 2.075, 2.05, 2.03, 2.02, 2.015, 1.985, 1.97, 1.955, 1.935,

		1.915, 1.91, 1.885, 1.86, 1.85, 1.835, 1.81, 1.795, 1.78, 1.765, 1.755, 1.73, 1.715, 1.69,

		1.675, 1.66, 1.645, 1.625, 1.61, 1.59, 1.57, 1.56, 1.54, 1.52, 1.5, 1.48, 1.46, 1.45, 1.43,

		1.41, 1.39, 1.37, 1.35, 1.33, 1.32, 1.3, 1.28, 1.26, 1.24, 1.22, 1.2, 1.18,

		1.16, 1.14, 1.12, 1.09, 1.07, 1.05, 1.03, 1.01, 0.985, 0.965, 0.94,

		0.91, 0.89, 0.86, 0.84, 0.8, 0.79, 0.76, 0.73, 0.7, 0.67, 0.64, 0.605,

		0.57, 0.535, 0.49, 0.45, 0.4, 0.345, 0.285, 0.2, 0.000 };

		int index;

		index = (int) (arg1 * 100);

		return (a[index]);

	}

	double bessel1(double arg) {

		// *** The solution of the Bessel-1 transcendent equation by dihotomia
		// ****

		double a, b, c;

		if (0 > arg)
			return 0;

		if ((Math.PI / 2) < arg) {
			return besj1(Math.PI / 2);
		}
		a = 0;

		b = 1.84;

		while (true) {
			c = (a + b) / 2;

			if (Math.abs(besj1(c) - arg) < 0.02) {
				return c;
			}

			if ((besj1(c) - arg) < 0) {
				a = c;
			} else {
				b = c;
			}

		}
	}

	double besj1(double arg) {

		// ***************** 1 order Bessel-function of the first kind
		// ***********

		double arg1 = 3 / arg, res;

		if (arg >= 3.0) {

			double f1, teta;

			teta = arg - 2.35619449 + 0.12499612 * (arg1) + 0.00005650 * Math.pow(arg1, 2) - 0.00637879
					* Math.pow(arg1, 3) + 0.00074348 * Math.pow(arg1, 4) + 0.00079824 * Math.pow(arg1, 5);

			f1 = 0.79788456 + 0.00000156 * (arg1) + 0.01659667 * Math.pow(arg1, 2) + 0.00017105 * Math.pow(arg1, 3)
					- 0.00249511 * Math.pow(arg1, 4) + 0.00113653 * Math.pow(arg1, 5);

			res = Math.pow(arg, -0.5) * f1 * Math.cos(teta);

			return res;

		}

		arg1 = arg / 3;

		res = 0.5 * arg - 0.56249985 * Math.pow(arg1, 2) * arg + 0.21093573 * Math.pow(arg1, 4) * arg - 0.03954289
				* Math.pow(arg1, 6) * arg + 0.00443319 * Math.pow(arg1, 8) * arg;

		return res;
	}

}
